Previous R Code | Home

# Setting up the environment, load the necessary packages
library(actuar)
library(dplyr)
library(EnvStats)
library(fitdistrplus)
library(kableExtra)
library(knitr)
library(ggplot2)
library(goftest)
library(MASS)
library(readxl)
library(tibble)
library(tidyr)

Loading Data

# load data
load(file = "census.RData")
load(file = "emissions.RData")
load(file = "hazards.RData")
load(file = "inflation.RData")
load(file = "rates.RData")
load(file = "prop_dist.RData")

# filtering non-impacting hazards
hazards <- hazards %>% arrange(Year) %>% filter(Property.Damage != 0)
# summarising hazard frequency
frequency <- hazards %>% group_by(Year) %>% summarise(events.pa = n())

Average Aggregate loss in n-years, using MLE Parameters

options(scipen=999)
set.seed(123)

n_years <- 10 # feel free to adjust: e.g. if you want aggregate loss in 5 years, put n_years <- 5
n_iterations <- 100000 # number of Monte_Carlo simulations: also feel free to adjust.
size <- 1.948557
mu <- 49.704918
shapelog <- 14.829296
ratelog <- 1.413602

sev <- c()

for (i in 1:n_iterations) {
  
  for (i in 1:n_years) {
    freq_vec <- rnbinom(n_years, size = size, mu = mu)
    sev_obs <- sum(rlgamma(freq_vec[i], shapelog = shapelog, ratelog = ratelog))
  }
  
  sev <- append(sev, sev_obs)
  
}

#sanity check: mean aggregate loss in n years / cumulative historical property damages: should be about = n_years / 60
(mean(sev)/sum(hazards$Property.Damage) )
## [1] 0.06110882
# Previusly we got 0.06579576. For n_years = 10, this should be about 0.1667

# Possible reason for small value: 
# the historical losses include some years where total losses are extreme, which isn't captured by just taking mean(sev)

sev_df <- data.frame(sev)
ggplot(sev_df, aes(x = sev)) + 
    geom_histogram(bins = 100) +
    xlim(c(0, 1000000000)) + ylim(c(0,15000)) +
    ggtitle("Aggregate Loss across Instances") +
    xlab("Aggregate Loss in N Years") +
    ylab("Sample")

Aggregate Loss per year, for each region

Prep

freq1 <- hazards %>% filter(Region == 1) %>% group_by(Year) %>% summarise(events.pa = n())
freq2 <- hazards %>% filter(Region == 2) %>% group_by(Year) %>% summarise(events.pa = n())
freq3 <- hazards %>% filter(Region == 3) %>% group_by(Year) %>% summarise(events.pa = n())
freq4 <- hazards %>% filter(Region == 4) %>% group_by(Year) %>% summarise(events.pa = n())
freq5 <- hazards %>% filter(Region == 5) %>% group_by(Year) %>% summarise(events.pa = n())
freq6 <- hazards %>% filter(Region == 6) %>% group_by(Year) %>% summarise(events.pa = n())

lsev1 <- hazards %>% filter(Region == 1) %>% mutate(Property.Damage = log(Property.Damage))
lsev2 <- hazards %>% filter(Region == 2) %>% mutate(Property.Damage = log(Property.Damage))
lsev3 <- hazards %>% filter(Region == 3) %>% mutate(Property.Damage = log(Property.Damage))
lsev4 <- hazards %>% filter(Region == 4) %>% mutate(Property.Damage = log(Property.Damage))
lsev5 <- hazards %>% filter(Region == 5) %>% mutate(Property.Damage = log(Property.Damage))
lsev6 <- hazards %>% filter(Region == 6) %>% mutate(Property.Damage = log(Property.Damage))

# fit historical frequency per region
nbin1 <- fitdist(freq1$events.pa, "nbinom")
nbin2 <- fitdist(freq2$events.pa, "nbinom")
nbin3 <- fitdist(freq3$events.pa, "nbinom")
nbin4 <- fitdist(freq4$events.pa, "nbinom")
nbin5 <- fitdist(freq5$events.pa, "nbinom")
nbin6 <- fitdist(freq6$events.pa, "nbinom")

# check fit using qq plots

all good!

# but frequency per year will increase because of climate change
# multiply the mu parameter by each RAF to get a table of different mus, depending on the climate scenario

freq_RAF1 <- emissions %>% mutate(
  SSP1..2.6 = SSP1..2.6*nbin1$estimate[[2]],
  SSP2.3.4 = SSP2.3.4*nbin1$estimate[[2]],
  SSP3.6.0 = SSP3.6.0*nbin1$estimate[[2]],
  SSP5.Baseline = SSP5.Baseline*nbin1$estimate[[2]],
  .keep = "unused"
)
freq_RAF2 <- emissions %>% mutate(
  SSP1..2.6 = SSP1..2.6*nbin2$estimate[[2]],
  SSP2.3.4 = SSP2.3.4*nbin2$estimate[[2]],
  SSP3.6.0 = SSP3.6.0*nbin2$estimate[[2]],
  SSP5.Baseline = SSP5.Baseline*nbin2$estimate[[2]],
  .keep = "unused"
)
freq_RAF3 <- emissions %>% mutate(
  SSP1..2.6 = SSP1..2.6*nbin3$estimate[[2]],
  SSP2.3.4 = SSP2.3.4*nbin3$estimate[[2]],
  SSP3.6.0 = SSP3.6.0*nbin3$estimate[[2]],
  SSP5.Baseline = SSP5.Baseline*nbin3$estimate[[2]],
  .keep = "unused"
)
freq_RAF4 <- emissions %>% mutate(
  SSP1..2.6 = SSP1..2.6*nbin4$estimate[[2]],
  SSP2.3.4 = SSP2.3.4*nbin4$estimate[[2]],
  SSP3.6.0 = SSP3.6.0*nbin4$estimate[[2]],
  SSP5.Baseline = SSP5.Baseline*nbin4$estimate[[2]],
  .keep = "unused"
)
freq_RAF5 <- emissions %>% mutate(
  SSP1..2.6 = SSP1..2.6*nbin5$estimate[[2]],
  SSP2.3.4 = SSP2.3.4*nbin5$estimate[[2]],
  SSP3.6.0 = SSP3.6.0*nbin5$estimate[[2]],
  SSP5.Baseline = SSP5.Baseline*nbin5$estimate[[2]],
  .keep = "unused"
)
freq_RAF6 <- emissions %>% mutate(
  SSP1..2.6 = SSP1..2.6*nbin6$estimate[[2]],
  SSP2.3.4 = SSP2.3.4*nbin6$estimate[[2]],
  SSP3.6.0 = SSP3.6.0*nbin6$estimate[[2]],
  SSP5.Baseline = SSP5.Baseline*nbin6$estimate[[2]],
  .keep = "unused"
)

# fit severity to log-gamma distribution

gamm1 <- fitdist(lsev1$Property.Damage, "gamma")
gamm2 <- fitdist(lsev2$Property.Damage, "gamma")
gamm3 <- fitdist(lsev3$Property.Damage, "gamma")
gamm4 <- fitdist(lsev4$Property.Damage, "gamma")
gamm5 <- fitdist(lsev5$Property.Damage, "gamma")
gamm6 <- fitdist(lsev6$Property.Damage, "gamma")

Simulation

n_iterations <- 10000 #number of Monte_Carlo simulations: 

# create empty vector for each region
sev1 <- c()
sev2 <- c()
sev3 <- c()
sev4 <- c()
sev5 <- c()
sev6 <- c()

for (i in 1:n_iterations) {
  freq_vec <- rnbinom(1, size = nbin1$estimate[[1]], mu = nbin1$estimate[[2]])
  sev_pa <- sum(exp(rgamma(freq_vec, shape = gamm1$estimate[[1]], rate = gamm1$estimate[[2]])))
  sev1 <- append(sev1, sev_pa)
}

for (i in 1:n_iterations) {
  freq_vec <- rnbinom(1, size = nbin2$estimate[[1]], mu = nbin2$estimate[[2]])
  sev_pa <- sum(exp(rgamma(freq_vec, shape = gamm2$estimate[[1]], rate = gamm2$estimate[[2]])))
  sev2 <- append(sev2, sev_pa)
}

for (i in 1:n_iterations) {
  freq_vec <- rnbinom(1, size = nbin3$estimate[[1]], mu = nbin3$estimate[[2]])
  sev_pa <- sum(exp(rgamma(freq_vec, shape = gamm3$estimate[[1]], rate = gamm3$estimate[[2]])))
  sev3 <- append(sev3, sev_pa)
}

for (i in 1:n_iterations) {
  freq_vec <- rnbinom(1, size = nbin4$estimate[[1]], mu = nbin4$estimate[[2]])
  sev_pa <- sum(exp(rgamma(freq_vec, shape = gamm4$estimate[[1]], rate = gamm4$estimate[[2]])))
  sev4 <- append(sev4, sev_pa)
}

for (i in 1:n_iterations) {
  freq_vec <- rnbinom(1, size = nbin5$estimate[[1]], mu = nbin5$estimate[[2]])
  sev_pa <- sum(exp(rgamma(freq_vec, shape = gamm5$estimate[[1]], rate = gamm5$estimate[[2]])))
  sev5 <- append(sev5, sev_pa)
}

for (i in 1:n_iterations) {
  freq_vec <- rnbinom(1, size = nbin6$estimate[[1]], mu = nbin6$estimate[[2]])
  sev_pa <- sum(exp(rgamma(freq_vec, shape = gamm6$estimate[[1]], rate = gamm6$estimate[[2]])))
  sev6 <- append(sev6, sev_pa)
}

# convert the filled in vector through iterations to dataframe
sev1_df <- data.frame(sev1)
sev2_df <- data.frame(sev2)
sev3_df <- data.frame(sev3)
sev4_df <- data.frame(sev4)
sev5_df <- data.frame(sev5)
sev6_df <- data.frame(sev6)

# such that we have the aggregate loss of each region for the following year as below

Value at Risk

#Percentiles: help compute VaR
perc <- c(0.25, 0.5, 0.75, 0.8, 0.85, 0.9, 0.95, 0.98, 0.99, 0.995, 0.999)
Percentile <- c(perc, "Mean", "SD", "Max")
gross1 <- c(quantile(sev1, probs = perc), mean(sev1), sd(sev1), max(sev1))
gross2 <- c(quantile(sev2, probs = perc), mean(sev2), sd(sev2), max(sev2))
gross3 <- c(quantile(sev3, probs = perc), mean(sev3), sd(sev3), max(sev3))
gross4 <- c(quantile(sev4, probs = perc), mean(sev4), sd(sev4), max(sev4))
gross5 <- c(quantile(sev5, probs = perc), mean(sev5), sd(sev5), max(sev5))
gross6 <- c(quantile(sev6, probs = perc), mean(sev6), sd(sev6), max(sev6))

results.table <- data.frame(Percentile, gross1, gross2, gross3, gross4, gross5, gross6)
# gross stands for the possible loss at the respective percentile
kable(results.table) %>%
  kable_styling(full_width = F)
Percentile gross1 gross2 gross3 gross4 gross5 gross6
0.25 771576.3 718714.2 443846.3 269850.7 141469.1 72640.08
0.5 3987297.7 4856075.9 3426386.3 1436337.3 948797.0 561336.72
0.75 18257501.2 22071997.9 18877478.1 6551006.8 4679887.4 3792230.23
0.8 26925007.0 30607123.6 28112198.9 9401245.0 6739647.3 6190111.53
0.85 45142458.0 45622136.8 45743680.9 14713714.1 10858227.3 11298922.94
0.9 86112377.5 77176865.7 84903025.7 25475294.9 19459361.5 23162083.99
0.95 242350621.9 181775554.4 227029497.3 60111612.4 54253928.1 69771423.59
0.98 852990634.1 504828181.8 704316318.7 216839921.7 193619072.3 311120070.60
0.99 2042810730.9 992463035.0 1596079613.4 512834674.1 511077369.7 969534745.57
0.995 5385551358.0 2413655738.3 5029179009.2 1023055975.5 1285980236.1 2314848111.77
0.999 38322605077.1 13295819215.7 22734842542.4 8286855343.2 7329127943.9 16392062142.30
Mean 5044036197.1 142664288.5 231297676.3 44281321.4 45027481.1 91931763.45
SD 482733286024.4 3902624225.4 8431638476.9 693474831.4 781354417.5 1933813793.56
Max 48271245840794.7 336517588734.2 801944277209.3 36685482548.7 42054941577.4 121258034304.42
# if we look at the historical property damage

hazards_sorted <- hazards %>% group_by(Year) %>% summarise(sum(Property.Damage))
head(hazards_sorted, 10)
## # A tibble: 10 × 2
##     Year `sum(Property.Damage)`
##    <int>                  <dbl>
##  1  1960                 303291
##  2  1961                 707416
##  3  1962                 238398
##  4  1963                1795433
##  5  1964                4623976
##  6  1965                2272604
##  7  1966                2473392
##  8  1967                2461031
##  9  1968                2225197
## 10  1969                1920886
# seems to underestimate aggregate loss
# may need to either exclude older years, or adjust damages for inflation

Incorporating Inflation

# formating the data for use
inflation = as.data.frame(inflation) %>% 
  remove_rownames %>% 
  column_to_rownames(var="Year")

# Data of inflation starts from 1962, we omit data before year 1962
hazards_adj <- hazards %>%
  dplyr::filter(Year >= 1962)

# reiterating the steps above with inflation incorporated
for (i in 1:nrow(hazards_adj)) {
  hazards_adj$Property.Damage[i] = hazards_adj$Property.Damage[i] * inflation$FV[rownames(inflation) == as.character(hazards_adj$Year[i])]
}

hazards <- hazards_adj

freq1 <- hazards %>% filter(Region == 1) %>% group_by(Year) %>% summarise(events.pa = n())
freq2 <- hazards %>% filter(Region == 2) %>% group_by(Year) %>% summarise(events.pa = n())
freq3 <- hazards %>% filter(Region == 3) %>% group_by(Year) %>% summarise(events.pa = n())
freq4 <- hazards %>% filter(Region == 4) %>% group_by(Year) %>% summarise(events.pa = n())
freq5 <- hazards %>% filter(Region == 5) %>% group_by(Year) %>% summarise(events.pa = n())
freq6 <- hazards %>% filter(Region == 6) %>% group_by(Year) %>% summarise(events.pa = n())

lsev1 <- hazards %>% filter(Region == 1) %>% mutate(Property.Damage = log(Property.Damage))
lsev2 <- hazards %>% filter(Region == 2) %>% mutate(Property.Damage = log(Property.Damage))
lsev3 <- hazards %>% filter(Region == 3) %>% mutate(Property.Damage = log(Property.Damage))
lsev4 <- hazards %>% filter(Region == 4) %>% mutate(Property.Damage = log(Property.Damage))
lsev5 <- hazards %>% filter(Region == 5) %>% mutate(Property.Damage = log(Property.Damage))
lsev6 <- hazards %>% filter(Region == 6) %>% mutate(Property.Damage = log(Property.Damage))

# fit historical frequency per region
nbin1 <- fitdist(freq1$events.pa, "nbinom")
nbin2 <- fitdist(freq2$events.pa, "nbinom")
nbin3 <- fitdist(freq3$events.pa, "nbinom")
nbin4 <- fitdist(freq4$events.pa, "nbinom")
nbin5 <- fitdist(freq5$events.pa, "nbinom")
nbin6 <- fitdist(freq6$events.pa, "nbinom")

# check fit using qq plots

looks all good!

# multiply the mu parameter by each RAF to get a table of different mus, depending on the climate scenario

freq_RAF1 <- emissions %>% mutate(
  SSP1..2.6 = SSP1..2.6*nbin1$estimate[[2]],
  SSP2.3.4 = SSP2.3.4*nbin1$estimate[[2]],
  SSP3.6.0 = SSP3.6.0*nbin1$estimate[[2]],
  SSP5.Baseline = SSP5.Baseline*nbin1$estimate[[2]],
  .keep = "unused"
)
freq_RAF2 <- emissions %>% mutate(
  SSP1..2.6 = SSP1..2.6*nbin2$estimate[[2]],
  SSP2.3.4 = SSP2.3.4*nbin2$estimate[[2]],
  SSP3.6.0 = SSP3.6.0*nbin2$estimate[[2]],
  SSP5.Baseline = SSP5.Baseline*nbin2$estimate[[2]],
  .keep = "unused"
)
freq_RAF3 <- emissions %>% mutate(
  SSP1..2.6 = SSP1..2.6*nbin3$estimate[[2]],
  SSP2.3.4 = SSP2.3.4*nbin3$estimate[[2]],
  SSP3.6.0 = SSP3.6.0*nbin3$estimate[[2]],
  SSP5.Baseline = SSP5.Baseline*nbin3$estimate[[2]],
  .keep = "unused"
)
freq_RAF4 <- emissions %>% mutate(
  SSP1..2.6 = SSP1..2.6*nbin4$estimate[[2]],
  SSP2.3.4 = SSP2.3.4*nbin4$estimate[[2]],
  SSP3.6.0 = SSP3.6.0*nbin4$estimate[[2]],
  SSP5.Baseline = SSP5.Baseline*nbin4$estimate[[2]],
  .keep = "unused"
)
freq_RAF5 <- emissions %>% mutate(
  SSP1..2.6 = SSP1..2.6*nbin5$estimate[[2]],
  SSP2.3.4 = SSP2.3.4*nbin5$estimate[[2]],
  SSP3.6.0 = SSP3.6.0*nbin5$estimate[[2]],
  SSP5.Baseline = SSP5.Baseline*nbin5$estimate[[2]],
  .keep = "unused"
)
freq_RAF6 <- emissions %>% mutate(
  SSP1..2.6 = SSP1..2.6*nbin6$estimate[[2]],
  SSP2.3.4 = SSP2.3.4*nbin6$estimate[[2]],
  SSP3.6.0 = SSP3.6.0*nbin6$estimate[[2]],
  SSP5.Baseline = SSP5.Baseline*nbin6$estimate[[2]],
  .keep = "unused"
)

# similarly also fit severity to log-gamma distribution

gamm1 <- fitdist(lsev1$Property.Damage, "gamma")
gamm2 <- fitdist(lsev2$Property.Damage, "gamma")
gamm3 <- fitdist(lsev3$Property.Damage, "gamma")
gamm4 <- fitdist(lsev4$Property.Damage, "gamma")
gamm5 <- fitdist(lsev5$Property.Damage, "gamma")
gamm6 <- fitdist(lsev6$Property.Damage, "gamma")

Inflation adjusted simulation

n_iterations <- 10000 #number of Monte_Carlo simulations: 

sev1 <- c()
sev2 <- c()
sev3 <- c()
sev4 <- c()
sev5 <- c()
sev6 <- c()
sev_agg <- numeric(n_iterations)

for (i in 1:n_iterations) {
  freq_vec <- rnbinom(1, size = nbin1$estimate[[1]], mu = nbin1$estimate[[2]])
  sev_pa <- sum(exp(rgamma(freq_vec, shape = gamm1$estimate[[1]], rate = gamm1$estimate[[2]])))
  sev1 <- append(sev1, sev_pa)
  sev_agg[i] = sev_agg[i] + sev_pa
}

for (i in 1:n_iterations) {
  freq_vec <- rnbinom(1, size = nbin2$estimate[[1]], mu = nbin2$estimate[[2]])
  sev_pa <- sum(exp(rgamma(freq_vec, shape = gamm2$estimate[[1]], rate = gamm2$estimate[[2]])))
  sev2 <- append(sev2, sev_pa)
  sev_agg[i] = sev_agg[i] + sev_pa
}

for (i in 1:n_iterations) {
  freq_vec <- rnbinom(1, size = nbin3$estimate[[1]], mu = nbin3$estimate[[2]])
  sev_pa <- sum(exp(rgamma(freq_vec, shape = gamm3$estimate[[1]], rate = gamm3$estimate[[2]])))
  sev3 <- append(sev3, sev_pa)
  sev_agg[i] = sev_agg[i] + sev_pa
}

for (i in 1:n_iterations) {
  freq_vec <- rnbinom(1, size = nbin4$estimate[[1]], mu = nbin4$estimate[[2]])
  sev_pa <- sum(exp(rgamma(freq_vec, shape = gamm4$estimate[[1]], rate = gamm4$estimate[[2]])))
  sev4 <- append(sev4, sev_pa)
  sev_agg[i] = sev_agg[i] + sev_pa
}

for (i in 1:n_iterations) {
  freq_vec <- rnbinom(1, size = nbin5$estimate[[1]], mu = nbin5$estimate[[2]])
  sev_pa <- sum(exp(rgamma(freq_vec, shape = gamm5$estimate[[1]], rate = gamm5$estimate[[2]])))
  sev5 <- append(sev5, sev_pa)
  sev_agg[i] = sev_agg[i] + sev_pa
}

for (i in 1:n_iterations) {
  freq_vec <- rnbinom(1, size = nbin6$estimate[[1]], mu = nbin6$estimate[[2]])
  sev_pa <- sum(exp(rgamma(freq_vec, shape = gamm6$estimate[[1]], rate = gamm6$estimate[[2]])))
  sev6 <- append(sev6, sev_pa)
  sev_agg[i] = sev_agg[i] + sev_pa
}
  
sev1_df <- data.frame(sev1)
sev2_df <- data.frame(sev2)
sev3_df <- data.frame(sev3)
sev4_df <- data.frame(sev4)
sev5_df <- data.frame(sev5)
sev6_df <- data.frame(sev6)
agg_df <- data.frame(sev_agg)

# such that we have the aggregate loss of each region for the following year as below

Value at Risk

#Percentiles: help compute VaR
perc <- c(0.25, 0.5, 0.75, 0.8, 0.85, 0.9, 0.95, 0.98, 0.99, 0.995, 0.999)
Percentile <- c(perc, "Mean", "SD", "Max")
gross1 <- c(quantile(sev1, probs = perc), mean(sev1), sd(sev1), max(sev1))
gross2 <- c(quantile(sev2, probs = perc), mean(sev2), sd(sev2), max(sev2))
gross3 <- c(quantile(sev3, probs = perc), mean(sev3), sd(sev3), max(sev3))
gross4 <- c(quantile(sev4, probs = perc), mean(sev4), sd(sev4), max(sev4))
gross5 <- c(quantile(sev5, probs = perc), mean(sev5), sd(sev5), max(sev5))
gross6 <- c(quantile(sev6, probs = perc), mean(sev6), sd(sev6), max(sev6))

results.table <- data.frame(Percentile, gross1, gross2, gross3, gross4, gross5, gross6)
# gross stands for the possible loss at the respective percentile
kable(results.table) %>%
  kable_styling(full_width = F)
Percentile gross1 gross2 gross3 gross4 gross5 gross6
0.25 1629122 1503668 695254.4 550895.2 333468.1 194038.8
0.5 7600693 8136864 5740506.6 2889949.3 1945553.3 1267204.2
0.75 31813657 33134941 31357535.1 11845188.0 9045598.4 6506298.2
0.8 45956634 45097667 46461285.3 16682941.2 12724498.7 9839231.7
0.85 69700881 66018839 73337046.9 25505572.0 19335621.8 15808223.5
0.9 118107966 107021259 132789273.4 42391769.1 35546309.5 31989814.9
0.95 299098477 238394085 314568412.3 92490508.8 85565374.3 84745089.9
0.98 969317650 638104307 960090973.5 256521679.0 258138884.0 327637895.6
0.99 2371394939 1247288513 2349725311.5 524540397.6 495050473.1 788602780.7
0.995 5249535516 2209462685 4772229097.5 1228020776.6 757472154.1 2059710239.8
0.999 40054422890 20203226080 32832237724.3 9355287842.4 2735246621.5 16217121377.4
Mean 1203891656 258056118 634654833.7 60393136.2 32302208.9 91630071.6
SD 99217385328 8705770893 30530976237.4 1149636007.6 375785525.7 1956144024.0
Max 9915532671009 704147416551 2766200940720.2 77977112578.9 26162819845.3 145163464873.0
# if we look at the historical property damage

hazards_sorted <- hazards %>% group_by(Year) %>% summarise(sum(Property.Damage))
head(hazards_sorted, 10)
## # A tibble: 10 × 2
##     Year `sum(Property.Damage)`
##    <int>                  <dbl>
##  1  1962               1999380.
##  2  1963              14743771.
##  3  1964              37223070.
##  4  1965              17926990.
##  5  1966              19027567.
##  6  1967              18327662.
##  7  1968              16129425.
##  8  1969              13455368.
##  9  1970              20713614.
## 10  1971              31066560.
# comparing the empirical and projected aggregate loss per region
emp_df1 <- data.frame(group = "emp", value = lsev1$Property.Damage)
emp_df2 <- data.frame(group = "emp", value = lsev2$Property.Damage)
emp_df3 <- data.frame(group = "emp", value = lsev3$Property.Damage)
emp_df4 <- data.frame(group = "emp", value = lsev4$Property.Damage)
emp_df5 <- data.frame(group = "emp", value = lsev5$Property.Damage)
emp_df6 <- data.frame(group = "emp", value = lsev6$Property.Damage)
emp_df <- data.frame(group = "emp", value = log(hazards_sorted$`sum(Property.Damage)`))

proj_df1 <- data.frame(group = "proj", value = log(sev1_df$sev1))
proj_df2 <- data.frame(group = "proj", value = log(sev2_df$sev2))
proj_df3 <- data.frame(group = "proj", value = log(sev3_df$sev3))
proj_df4 <- data.frame(group = "proj", value = log(sev4_df$sev4))
proj_df5 <- data.frame(group = "proj", value = log(sev5_df$sev5))
proj_df6 <- data.frame(group = "proj", value = log(sev6_df$sev6))
proj_df <- data.frame(group = "proj", value = log(agg_df$sev_agg))

> seems to underestimate aggregate loss.
> may need to either exclude older years, or adjust damages for inflation.

Program Loss

Assumptions

# Set up parameters
n_iterations <- 10000 # number of Monte_Carlo simulations: 
set.seed(123)
mat_increase <- runif(n_iterations, 0,0.5)
hh_cost <- runif(n_iterations, 0.4, 0.75)

Short-term Program Cost Projection

# We omit region 5 and 6 as they are low-risk region and would not relocate under our program design
loss1 <- c()
loss2 <- c()
loss3 <- c()
loss4 <- c()
loss_agg <- numeric(n_iterations)

for (i in 1:n_iterations) {
  freq_vec <- rnbinom(1, size = nbin1$estimate[[1]], mu = nbin1$estimate[[2]])
  sev_pa <- sum(exp(rgamma(freq_vec, shape = gamm1$estimate[[1]], rate = gamm1$estimate[[2]])))
  loss <- sev_pa * (1 + mat_increase[i]) * (1 + hh_cost[i])* prop_dist[1,1]
  loss1 <- append(loss1,loss)
  loss_agg[i] = loss_agg[i] + loss
}

for (i in 1:n_iterations) {
  freq_vec <- rnbinom(1, size = nbin2$estimate[[1]], mu = nbin2$estimate[[2]])
  sev_pa <- sum(exp(rgamma(freq_vec, shape = gamm2$estimate[[1]], rate = gamm2$estimate[[2]])))
  loss <- sev_pa * (1 + mat_increase[i]) * (1 + hh_cost[i])* prop_dist[2,1]
  loss2 <- append(loss2,loss)
  loss_agg[i] = loss_agg[i] + loss
}

for (i in 1:n_iterations) {
  freq_vec <- rnbinom(1, size = nbin3$estimate[[1]], mu = nbin3$estimate[[2]])
  sev_pa <- sum(exp(rgamma(freq_vec, shape = gamm3$estimate[[1]], rate = gamm3$estimate[[2]])))
  loss <- sev_pa * (1 + mat_increase[i]) * (1 + hh_cost[i])* prop_dist[3,1]
  loss3 <- append(loss3,loss)
  loss_agg[i] = loss_agg[i] + loss
}

for (i in 1:n_iterations) {
  freq_vec <- rnbinom(1, size = nbin4$estimate[[1]], mu = nbin4$estimate[[2]])
  sev_pa <- sum(exp(rgamma(freq_vec, shape = gamm4$estimate[[1]], rate = gamm4$estimate[[2]])))
  loss <- sev_pa * (1 + mat_increase[i]) * (1 + hh_cost[i])* prop_dist[4,1]
  loss4 <- append(loss4,loss)
  loss_agg[i] = loss_agg[i] + loss
}

perc <- c(0.25, 0.5, 0.75, 0.8, 0.85, 0.9, 0.95, 0.98, 0.99, 0.995, 0.999)
Percentile <- c(perc, "Mean", "SD", "Max")
gross1 <- c(quantile(loss1, probs = perc), mean(loss1), sd(loss1), max(loss1))
gross2 <- c(quantile(loss2, probs = perc), mean(loss2), sd(loss2), max(loss2))
gross3 <- c(quantile(loss3, probs = perc), mean(loss3), sd(loss3), max(loss3))
gross4 <- c(quantile(loss4, probs = perc), mean(loss4), sd(loss4), max(loss4))

results.table2 <- data.frame(Percentile, gross1, gross2, gross3, gross4)
kable(results.table2) %>%
  kable_styling(full_width = F)
Percentile gross1 gross2 gross3 gross4
0.25 1773280 1019100 74529.72 83517.7
0.5 8360132 5483970 600852.80 428597.6
0.75 32746798 22139960 3137050.00 1734636.5
0.8 46909268 30821064 4602860.87 2439286.1
0.85 70005650 45662599 6952049.40 3690976.2
0.9 125538099 77121843 12422503.49 6287813.5
0.95 302489888 172377322 30382419.13 13516231.4
0.98 843985864 480661320 76733338.24 39050930.3
0.99 2003026188 1023430145 158750236.36 74531533.8
0.995 4972522095 2013416270 287780740.88 138234940.7
0.999 27674584845 11762993830 1151862498.42 867045607.5
Mean 558710983 95374114 14223755.12 6267264.7
SD 39654558069 1467784285 278623079.60 76366550.0
Max 3960776747798 69502120999 23841760582.85 4617336066.4

Long-term Loss Projection (2030)

sev1 <- c()
sev2 <- c()
sev3 <- c()
sev4 <- c()
sev5 <- c()
sev6 <- c()
sev_agg <- numeric(n_iterations)


for (i in 1:n_iterations) {
  freq_vec <- rnbinom(1, size = nbin1$estimate[[1]], mu = freq_RAF1[2,4])
  sev_pa <- sum(exp(rgamma(freq_vec, shape = gamm1$estimate[[1]], rate = gamm1$estimate[[2]])))
  sev1 <- append(sev1, sev_pa)
  sev_agg[i] = sev_agg[i] + sev_pa
}

for (i in 1:n_iterations) {
  freq_vec <- rnbinom(1, size = nbin2$estimate[[1]], mu = freq_RAF2[2,4])
  sev_pa <- sum(exp(rgamma(freq_vec, shape = gamm2$estimate[[1]], rate = gamm2$estimate[[2]])))
  sev2 <- append(sev2, sev_pa)
  sev_agg[i] = sev_agg[i] + sev_pa
}

for (i in 1:n_iterations) {
  freq_vec <- rnbinom(1, size = nbin3$estimate[[1]], mu = freq_RAF3[2,4])
  sev_pa <- sum(exp(rgamma(freq_vec, shape = gamm3$estimate[[1]], rate = gamm3$estimate[[2]])))
  sev3 <- append(sev3, sev_pa)
  sev_agg[i] = sev_agg[i] + sev_pa
}

for (i in 1:n_iterations) {
  freq_vec <- rnbinom(1, size = nbin4$estimate[[1]], mu = freq_RAF4[2,4])
  sev_pa <- sum(exp(rgamma(freq_vec, shape = gamm4$estimate[[1]], rate = gamm4$estimate[[2]])))
  sev4 <- append(sev4, sev_pa)
  sev_agg[i] = sev_agg[i] + sev_pa
}

for (i in 1:n_iterations) {
  freq_vec <- rnbinom(1, size = nbin5$estimate[[1]], mu = freq_RAF5[2,4])
  sev_pa <- sum(exp(rgamma(freq_vec, shape = gamm5$estimate[[1]], rate = gamm5$estimate[[2]])))
  sev5 <- append(sev5, sev_pa)
  sev_agg[i] = sev_agg[i] + sev_pa
}

for (i in 1:n_iterations) {
  freq_vec <- rnbinom(1, size = nbin6$estimate[[1]], mu = freq_RAF6[2,4])
  sev_pa <- sum(exp(rgamma(freq_vec, shape = gamm6$estimate[[1]], rate = gamm6$estimate[[2]])))
  sev6 <- append(sev6, sev_pa)
  sev_agg[i] = sev_agg[i] + sev_pa
}

sev1_df <- data.frame(sev1)
sev2_df <- data.frame(sev2)
sev3_df <- data.frame(sev3)
sev4_df <- data.frame(sev4)
sev5_df <- data.frame(sev5)
sev6_df <- data.frame(sev6)
agg_df <- data.frame(sev_agg)

perc <- c(0.25, 0.5, 0.75, 0.8, 0.85, 0.9, 0.95, 0.98, 0.99, 0.995, 0.999)
Percentile <- c(perc, "Mean", "SD", "Max")
gross1 <- c(quantile(sev1, probs = perc), mean(sev1), sd(sev1), max(sev1))
gross2 <- c(quantile(sev2, probs = perc), mean(sev2), sd(sev2), max(sev2))
gross3 <- c(quantile(sev3, probs = perc), mean(sev3), sd(sev3), max(sev3))
gross4 <- c(quantile(sev4, probs = perc), mean(sev4), sd(sev4), max(sev4))
gross5 <- c(quantile(sev5, probs = perc), mean(sev5), sd(sev5), max(sev5))
gross6 <- c(quantile(sev6, probs = perc), mean(sev6), sd(sev6), max(sev6))

results.table3 <- data.frame(Percentile, gross1, gross2, gross3, gross4, gross5, gross6)
kable(results.table3) %>%
  kable_styling(full_width = F)
Percentile gross1 gross2 gross3 gross4 gross5 gross6
0.25 2253827 2027259 1082091 841328.6 475780.7 283785.4
0.5 9871308 10154043 8430743 3914267.2 2599826.3 1711185.4
0.75 39292323 40692292 41054365 15310939.5 10362078.5 8688070.7
0.8 55428815 56744851 58447877 20933642.9 14602333.8 13313130.0
0.85 81825775 81338836 92664312 31407893.0 22378339.4 22406803.2
0.9 140583587 132491259 166252708 54261843.8 38806285.9 41667822.9
0.95 324909592 277850366 399773690 120193097.4 92573237.3 119270332.0
0.98 1062707651 708535351 1281936291 357693151.0 251844867.1 438675650.5
0.99 2514425435 1554361731 2636870598 686498785.9 560331425.0 1147060237.3
0.995 4907374494 3551831203 4658265143 1535520349.1 1434040640.1 2740432992.3
0.999 26336456382 13755101965 19262003665 7942801272.7 10588742218.8 15540346011.2
Mean 267774620 183964977 322341756 105824411.7 149074479.4 1361449511.9
SD 8503374319 5102168758 16347868942 4541330211.3 8513156508.6 123506949096.5
Max 778764582297 475226499381 1624585172773 428792212121.5 843691653666.6 12344106818691.6

Long-term Program Cost Projection

mat_increase <- runif(n_iterations, 0,0.5)
hh_cost <- runif(n_iterations, 0.4, 0.75)

loss1 <- c()
loss2 <- c()
loss3 <- c()
loss4 <- c()
loss_agg <- numeric(n_iterations)

for (i in 1:n_iterations) {
  freq_vec <- rnbinom(1, size = nbin1$estimate[[1]], mu = freq_RAF1[2,4])
  sev_pa <- sum(exp(rgamma(freq_vec, shape = gamm1$estimate[[1]], rate = gamm1$estimate[[2]])))
  loss <- sev_pa * (1 + mat_increase[i]) * (1 + hh_cost[i])* prop_dist[1,1]
  loss1 <- append(loss1,loss)
  loss_agg[i] = loss_agg[i] + loss
}

for (i in 1:n_iterations) {
  freq_vec <- rnbinom(1, size = nbin2$estimate[[1]], mu = freq_RAF2[2,4])
  sev_pa <- sum(exp(rgamma(freq_vec, shape = gamm2$estimate[[1]], rate = gamm2$estimate[[2]])))
  loss <- sev_pa * (1 + mat_increase[i]) * (1 + hh_cost[i])* prop_dist[2,1]
  loss2 <- append(loss2,loss)
  loss_agg[i] = loss_agg[i] + loss
}

for (i in 1:n_iterations) {
  freq_vec <- rnbinom(1, size = nbin3$estimate[[1]], mu = freq_RAF3[2,4])
  sev_pa <- sum(exp(rgamma(freq_vec, shape = gamm3$estimate[[1]], rate = gamm3$estimate[[2]])))
  loss <- sev_pa * (1 + mat_increase[i]) * (1 + hh_cost[i])* prop_dist[3,1]
  loss3 <- append(loss3,loss)
  loss_agg[i] = loss_agg[i] + loss
}

for (i in 1:n_iterations) {
  freq_vec <- rnbinom(1, size = nbin4$estimate[[1]], mu = freq_RAF4[2,4])
  sev_pa <- sum(exp(rgamma(freq_vec, shape = gamm4$estimate[[1]], rate = gamm4$estimate[[2]])))
  loss <- sev_pa * (1 + mat_increase[i]) * (1 + hh_cost[i])* prop_dist[4,1]
  loss4 <- append(loss4,loss)
  loss_agg[i] = loss_agg[i] + loss
}

perc <- c(0.25, 0.5, 0.75, 0.8, 0.85, 0.9, 0.95, 0.98, 0.99, 0.995, 0.999)
Percentile <- c(perc, "Mean", "SD", "Max")
gross1 <- c(quantile(loss1, probs = perc), mean(loss1), sd(loss1), max(loss1))
gross2 <- c(quantile(loss2, probs = perc), mean(loss2), sd(loss2), max(loss2))
gross3 <- c(quantile(loss3, probs = perc), mean(loss3), sd(loss3), max(loss3))
gross4 <- c(quantile(loss4, probs = perc), mean(loss4), sd(loss4), max(loss4))

results.table4 <- data.frame(Percentile, gross1, gross2, gross3, gross4)
kable(results.table4) %>%
  kable_styling(full_width = F)
Percentile gross1 gross2 gross3 gross4
0.25 2335259 1342749 100769.9 112353.7
0.5 10355933 7051089 758106.0 541429.2
0.75 41302714 26843289 3927277.9 2203397.0
0.8 57970835 37150681 5635489.4 3011946.9
0.85 88649842 53646685 8873599.8 4615253.1
0.9 151413958 87399398 15567404.2 7691604.5
0.95 367651779 187222273 39951542.0 17245138.2
0.98 1161887954 474111273 116524114.7 49251027.9
0.99 2547962197 922411230 272785822.1 106019873.7
0.995 5908277721 2062009198 624060387.7 205903684.5
0.999 37035206478 5409505917 2915694448.3 824002405.7
Mean 424106761 63537714 21106645.3 9693366.6
SD 19462748071 433677206 279353331.3 218267614.1
Max 1917033156560 22233498667 17929860326.9 19770385742.4

Previous R Code | Home